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Abstract 


A computer  program  aimed  at  the  phase  separation  between  gas  and  liquid  at 
zero  gravity,  induced  by  vortex  motion,  is  developed.  It  utilizes  an  explicit 
solution  method  for  a set  of  equations  describing  rotating  gas-liquid  flows. 

The  vortex  motion  is  established  by  a tangential  fluid  injection.  A Lax- 
Wendroff  two-step  (McCormack's)  numerical  scheme  is  used.  This  program  can  be 
used  to  study  the  fluid  dynamical  behavior  of  the  rotational  two-phase  fluids  in 
a cylindrical  tank.  It  provides  a quick/easy  sensitivity  test  on  various 
parameters  and  thus  provides  the  guidance  for  the  design  and  use  of  actual 
physical  systems  for  handling  two-phase  fluids. 

Key  Words:  computer  code;  gas-liquid  separation;  numerical  modeling;  two-phase 

vortex  motions 
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I.  Introduction 


Mechanical  systems  have  been  devised  for  producing  artificial  gravity 
fields  to  spin-up  liquids  in  containers.  These  involve  rotating  mechanisms  that 
are  cumbersome  and,  more  importantly,  have  moving  parts  that  can  wear  out. 

Here,  liquid  rotation  created  by  fluid  injection  is  considered.  The  detail 
analysis  of  the  two-phase  vortex  model  can  be  found  elsewhere  [1].  In  this 
report,  the  details  of  the  computer  code  are  described. 

The  computer  program  was  developed  to  study  the  fluid  dynamical  behavior  of 
two-phase  fluids  in  a tank  at  zero  gravity.  The  phase  separation  between  gas 
and  liquid,  induced  by  vortex  motions,  is  of  primary  interest.  The  program 
utilizes  an  explicit  solution  method  for  a set  of  equations  describing  rotating 
gas-liquid  flows.  The  vortex  motion  is  established  by  a tangential  fluid 
injection.  A Lax-Wendroff  two-step  (McCormack's)  numerical  scheme  is  used  in 
the  computer  program.  This  scheme  uses  a conservation  form  of  a system  of 
equations  together  with  an  auto  time  step  feature. 

The  program  was  developed  and  tested  on  an  HP-1000  minicomputer.  The  HP- 
1000  's  FORTRAN  77  is  based  on  the  American  National  Standards  Institute  (ANSI) 

77  standard  programming  language  FORTRAN  (ANSI  X3. 9-1978).  The  HP  FORTRAN  77 
has  extensions  to  provide  a more  structured  approach  to  program  development  and 
more  flexibility  in  computing  for  scientific  applications.  It  fully  implements 
the  Military  Standard  Definition  (MIL-STD-1 753)  of  extensions  to  the  ANSI  77 
standard.  In  order  to  make  the  computer  code  more  useful  for  other  computer 
systems,  modifications  have  been  made  so  that  the  code  is  closer  to  the  ANSI  77 
standard  and  thus  less  system  dependent.  Some  limited  extensions  are  still  kept 
in  order  to  produce  the  code  in  the  HP-1000.  Since  the  graphic  routines  are 
system  dependent  and  must  be  modified  with  their  equivalents  at  each  computing 
facility,  the  original  graphic  code  has  not  been  included  in  this  report.  All 


lines  preceded  by  ”*V'’  are  originally  adopted  to  use  the  vector  operation 
package  supplied  by  Hewlett  Packard.  The  speed  of  the  code  can  be  increased  by 
replacing  many  ”do-loop"  operations  in  the  code  with  high  speed  vector 
operations.  Little  effort  is  required  to  incorporate  the  vector  operation  into 
the  code  if  the  vector  operation  package  is  in  the  system. 

Thus  with  limited  effort,  the  program  can  be  adapted  easily  to  most  systems 
accepting  the  ANSI  77  standard  FORTRAN.  For  example,  the  EMA  (Extended  Memory 
Area)  statements  may  have  to  be  removed  from  the  code  for  some  computers.  Also 
double  precision  real  numbers  (Real  * 8)  could  be  replaced  by  single  precision 
real  numbers. 


II . Model  Equations 

The  vortex  induced  model  is  based  on  a two-phase,  two-fluid  continuum  [2]. 
It  incorporates  several  interactions  between  phases;  namely  fluid  drag  and 
virtual  mass  effects  and  it  can  be  modified  to  include  additional  interaction 
effects.  Detailed  analysis  of  the  model  has  been  reported  in  Ref.  1.  A brief 
summary  of  the  system  of  equations  is  given  below. 

The  equations  for  the  conservation  of  mass  and  momentum  for  the  two  fluid 
two-phase  model  in  an  one-dimensional,  axisymmetric  case 
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^d2  = > 


and 


<p  > = a^a^p^p2  + A^  ^ <^2^2^ 


The  effective  stresses  are  modeled  as 


T , = 2y,  3V  , /3r 
rrk  k rk 


a,  ^ ■’^Q  1 =1^1  ^ /r)/3r 
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X = 2 u V /r 
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with 

e t 

^k  = "k  " \ 

and  the  interfacial  forces  are  modeled  in  the  form  of 

Ml  - A^(?2  - Vi)  * ^ (V^  - ?i) 


is  the  force  density  acting  on  the  phase  1 by  the  phase  2.  and  are  the 

added  mass  and  drag  coefficients,  respectively. 

The  incompressibility  condition  is  reduced  to  ■*■  ^2^r2  ” 

Q is  the  net  radial  outflow, 
r 

In  the  program  = 0 is  assumed,  since  the  mixture  pumped  out  Is  Injected 
immediately  back  into  the  tank  at  the  nearby  location.  The  net  volume  or  mass 
in  the  system  is  effectively  unchanged  except  for  the  net  change  on  the  angular 
momentum.  Thus,  the  pump  system  (withdrawal  and  injection)  acts  as  a body  for-'e 
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on  the  mixture  at  the  nozzle  location.  The  net  momentum  gain  is  thus  the 
momentum  introduced  into  the  system  minus  the  local  momentum  pumped  out. 
Therefore,  we  will  model  this  pumping  dynamic  by  body  forces  without  considering 
the  mass  transfer.  That  is,  the  body  force  density  will  be  replaced  by 


a.  p,  V . 

k j — — 

the  net  momentum  gain,  — r — - (V.n  - V) 

2ir  j 

injection  speed. 


at  the  nozzle  location. 


where  V . 

J 


is  the 


III . Numerical  Method 

The  complete  solution  of  the  complicated  system  of  equations  can  only  be 
obtained  through  numerical  methods.  An  improved  Lax-Wendrof f , two-step  scheme, 
(also  referred  to  as  MacCormack’s  method)  [3,  4]  is  adopted  for  solving  this 
time-dependent  problem.  This  non-centere<f  differencing  scheme,  using  a full 
step  backward  prediction  and  forward  correction  version,  requires  no  explicit 
artificial  viscosity  if  a proper  stability  condition  is  satisfied.  Using  this 
technique  for  solving  fluid  flow  problems  is  very  efficient  and  has  been  in 
widespread  and  successful  use  for  some  time.  It  is  good  both  for  the  time- 
accurate  computation  of  steady  and  unsteady  flow  problems.  The  general  features 
of  the  scheme  are:  i)  its  explicitly  conservative  form,  ii)  it  is  a two-step 

predictor-correction  type,  iii)  it  is  three  point,  two  level  - that  is,  the 

n+1  n 

solution  of  f^  at  level  n+1  depends  only  on  three  values  of  f^  at  level  n,  and 

iv)  it  is  second-order  accurate  in  time  and  in  space. 

For  using  the  MacCormack's  numerical  technique,  the  system  of  equations  can 

be  expressed  in  the  conservative  form  as: 

W = F + P + gG  + S 
t r r r 

Here  the  subscripts  (t  and  r)  denote  partial  differentiation  with  respect  to  t 
and  r,  respectively,  and  W,  F,  P^,  gG^  and  S are  column  matrices  with  five 
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elements.  All  the  components  of  F,  P , gG  and  S can  be  regarded  as  functions 

r r 

of  the  components  of  W which  are  the  independent  variables.  The  fundamental 
theory  of  the  MacCormack’s  scheme  is  briefly  given  below. 

For  second  order  accuracy,  the  solution  could  be  written  as 

= W®  + At  W ° W ° 


,,o  At  ,,  o ^ At  0 ^ o. 

= W + — W.  + :r-  (W.  + At  W.  . ) 

2 t 2 t tt 


where 


and 


= I (W°  + At  W^°)  + I (W°  + At 

= I (wP  + W°) 

+ At  is  the  predicted  value, 

= W°  + At  W ^ is  the  corrected  value. 


The  superscripts  denote  the  time-level  of  the  information  and  subscripts  denote 
the  partial  derivative  with  respect  to  either  time  t or  space  r.  Specifically, 
superscripts  0 and  1 are  the  initial  and  the  completely  advanced  time  (here  two 
steps)  plane;  p and  c are  the  predicted  (1st  step)  and  corrected  (2nd  step)  time 
plane.  Thus,  W ^ is  the  time  derivative  of  W evaluated  at  the  initial  time,  and 
is  time  derivative  of  W evaluated  at  the  predicted  time. 

Fig.  1 shows  the  diagram  of  the  two  step  difference  scheme  used  in  the 
computer  program.  Due  to  the  difference  scheme,  the  spatial  location  after  each 
step  in  time  is  a half  grid  off  from  the  original  one.  Thus,  the  spatial  offset 
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which  resulted  from  a backward  predicting  step  will  cancel  with  those  of  the 
forward  correcting  step. 

Numerically,  the  predicted  values  are 


= i.  ) 

”i  2 ^”i-l/2  '^i+1/2^ 


where 


= i (W.°^  . W.°)  . [(F.°  - F.°^)  . 


. At  (3O  , 3^0^,  , 


and  the  corrected  value  is  evaluated  at  the  predicted  time  place,  that  is  at 


“in/2- 


w“ 

L 


W.°  + At 
1 t 


= W ° + — [(F  ^ 

i Ar  i+1/2 


F ^ ) + 

i-1/2^ 


(g  ^ + g ^ ) 

^®i  + l/2  *^i-l/2^ 


<h*i/2  - 


+ — (S^  +3^^  )+AtP^ 

2 ^^i+1/2  ^i-1/2''  ^i 


P c 

Here  pressure  correction  terms  at  half  and  full  time 


steps  respectively.  Thus,  for  each  time  step,  the  advance  is  carried  out  in  two 
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steps:  a full  step  backward  predictor,  and  then  a forward  corrector.  As 

indicated  in  the  diagram,  the  subscript  i is  the  regular  mesh  spatial  location 
at  which  solution  is  to  be  advanced,  i _+  1 is  the  spatial  location  of  regular 
mesh  points  immediately  to  the  right  and  left  of  the  location  i,  i + 1/2  is  the 
location  midway  between  i and  i + 1 or  between  i - 1 and  i at  the  predictor 
plane.  Thus,  for  each  time  step  as  the  procedure  advanced,  the  outermost  data 
points  at  the  boundary  are  not  updated  through  the  numerical  scheme.  The  values 
at  the  boundary  are  to  be  given  through  some  suitable  boundary  conditions.  The 
numerical  procedure  utilizes  a uniformly  preselected  spatial  mesh  and  variable 
time  increment.  To  avoid  a singularity  at  the  center  of  the  core  region,  a 

finite  radius  is  used  for  the  inner  boundary.  The  tank  radius  R is  the  outer 

boundary.  The  time  step  is  determined  at  each  time  step  to  ensure  numerical 

stability  [5].  For  a finite  grid  size  Ar,  the  maximum  time  step  At  is  given  by 


At,  = 1/[|C 


dk 


rk 


/Ar  + 


Ar 


“2*'2 


=.2>’ 


Where  k = 1 and  2.  The  minimum  At^^  (with  some  rounding  off)  is  used  for  the 
time  step.  Normally,  the  technique  with  the  time  step  condition  gives  fairly 
good  numerical  stability.  However,  in  critical  conditions  numerical  damping  can 
be  added  either  for  damping  oscillations  due  to  large  gradients  or  for 
accelerating  the  calculation  by  increasing  the  time  step.  A damping  factor,  D 
thus  was  added  in  the  program  as 


W. 

1 


ID 


(1-D)  + (W^_^  + 


D 
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1_  XD 

where  is  the  value  obtained  based  on  the  two-step  scheme,  and  W is  the 

value  after  the  damping  factor  D is  added.  A typical  value  of  D = 0.2  can  be 

used  for  debugging  the  program.  If  no  damping  factor  is  desired,  D = 0 should 

be  used. 

The  computer  program  was  written  in  a Fortran  77  based  computer  code.  The 
code  will  permit  evaluation  of  the  effects  of  various  parameters  which  control 
the  fluid  dynamical  behavior.  These  include  tank  size,  fluid  properties,  such 
as  density  and  viscosity,  etc.,  characteristic  gas  bubble  and  liquid  drops 
sizes,  and  relative  location  of  injection  nozzles. 

A sample  input  and  its  output  are  shown  on  Exhibits  A and  3,  respectively. 
The  initial  conditions  for  the  gas  and  liquid  volume  fractions  are  taken  to  be 
25J  gas  and  75J  liquid.  These  fractions  are  uniformly  distributed  over  the 
circular  cross-section  of  the  cylindrical  tank.  Initially  both  fluids  are  at 
rest.  Other  parameters  can  be  found  in  the  sample  input  in  Exhibit  A.  The 
resulting  velocity  distributions  and  gas  volume  fraction  as  function  of  time  for 
the  sample  run  are  given  in  Figs.  (2)  and  (3)  respectively.  The  velocity 
distributions  are  displayed  along  equally  spaced  rays  at  different  times  to 
enable  clear  observation.  These  velocity  vector  fields  indicate  all  flows  are 
primarily  in  angular  rotation  with  gas  phase  tending  to  move  inward  and  liquid 
phase  trying  to  move  outward,  as  expected.  As  the  result  of  these  radial 
movements,  the  volume  fraction  distribution  is  also  changed  with  time.  And  as 
expected,  the  gas  volume  fraction  is  increasing  at  the  inner  region  and 
decreasing  at  the  outer  region  as  shown  in  Figure  3.  More  detailed  results  have 
been  reported  in  Ref.  1 . 


9 


IV.  Program  Details 


The  complete  computer  code  is  listed  in  Appendix  A.  The  code  consists  of  a 
main  program,  GLVM  and  several  subroutines.  It  is  written  in  subroutine  form 
such  that  each  subroutine  performs  an  individual  task.  Each  logical  part  is 
clearly  isolated  and  it  can  be  easily  modified  to  reflect  different  modelings 
for  the  interfacial  forces.  The  interactive  input  mode  with  self-instruction  is 
used  for  easy  parameter  insertion.  Many  instructive  internal  documentations  are 
included  in  the  program.  In  the  following,  each  subroutine  is  listed  with  a 
brief  description  of  its  major  function. 

1)  GLVM,  the  Main  Program. 

*To  initialize  data  and  start  the  program:  Logical  unit  to  save  data  (LUS), 
(Logical  Unit  is  1 for  terminal,  and  6 for  printer),  job  identification 
notes  (NOTES) i data  file  name  for  saving  data  (NAMR),  initial  time  (TO), 
final  time  (TMAX) , time  interval  for  data  output  (DTPRT),  etc. 

*To  control  the  calling  sequences  to  the  other  subroutines. 

*To  check  the  time  step. 

*To  save,  print  (and  plot)  the  output  data. 

*To  obtain  the  predicted  and  corrected  values  in  the  two  step,  numerical 
scheme. 

*To  impose  boundary  condition. 

*To  update  the  data,  time,  and  step  number  for  the  new  time  step. 

*To  provide  a shutdown  procedure  either  in  normal  (e.g.,  t > t ) or 

max 

abnormal  (e.g..  At  is  too  small)  conditions. 

2)  INIT,  Initialization. 

*To  input  the  test  parameters,  initial  conditions  and  set-up  the  initial 
column  matrix  W. 
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Default  values  are  provided  for  most  of  the  parameters.  The  default  values 
are  listed  at  each  interactive  input  step.  If  the  default  value  is  acceptable, 
a comma  is  inputted. 

Some  of  the  relevant  symbols  used  are  listed  below: 


ALMT 

Limit  values  of  a^,  ALMT(l)  < < ALMT(2). 

DAMP 

Numerical  damping  factor.  Normally  set  to  0. 

DEN1D2 

Density  ratio,  P^_/p2 

DO 

Y 

Base  diameter,  i.e.,  d,  = d , a, 

k ok  k 

DS 

Density  scale  = 

EVF 

effective  eddy  viscosity  factor. 

GAMMA 

Diameter  exponent  Y 

IVTX 

Type  of  simple  initial  flow:  0 = at  rest,  1 = simple  rotation,  2 = 
H.amrael-Oseen  Vortex,  3 = G.I.  Taylor  Vortex.  This  is  needed  only  when 
there  is  no  data  file  (NAMR)  given  for  an  initial  condition. 

IW 

Boundary  condition  at  wall  (for  tangential  component).  0 = free-slip 

(no-skin  friction),  1 = non-slip.  Also,  when  I = 1,  a factor  of  (1- 

w 

r)*^’^  was  included  on  IVTX  flow  to  simulate  an  initial  power  law 
boundary  layer. 

MM 

Size  of  data  arrays,  MM  > NG.  MM  appears  in  many  subroutines. 

MU 

dynamic  viscosity  for  phase  k. 

MUEF 

(1  + EVF)  * MU.  Effective  viscosity. 

MU1D2 

^l'^^2’  Viscosity  ratio. 

NA,  ND 

n ,n  , Exponents  for  weighting  function  for  drag  and  added  mass 

3.  D 

coefficients,  A and  A,. 

a d 

NAMR 

Data  file  name  for  initial  condition,  if  any. 

NG 

Number  of  grid  points  used. 

OMEGA 

0),  initial  rotation  speed,  V.  = or,  when  IVTX  > 0. 

y 

PS 

2 

Pressure  scale,  DS  * VS 

QJ,VJ 

Injection  flow  rate  and  speed. 

RE 

Reynolds  number,  p V R/u- 

w 3 ^ 

RJ1,RJ2 

Jet  opening,  RJl  < r < RJ2 

RPEAK 

Location  of  the  peak  speed  of  the  initial  vortex,  if  IVTX  > 1 

RTANK 

Tank  radius  = Length  scale  LS. 

VPEAK 

Peak  speed  of  the  initial  vortex,  if  IVTX  > 1. 

TJ1,TJ2 

Injection  time,  TJl  < t < TJ2. 

TS 

Time  scale  = LS/VS 

VS 

Velocity  scale. 

The 

format  of  the  data  file  for  the  initial  condition  (if  any)  is  a six 

column  and  NG  row  data  file,  where  NG  is  the  number  of  grid  points.  The  column 
sequence  is  K,  a^(K),  ^91^^^’  '^92^^^'  where  K is  the  grid  point 

number,  and  the  rest  of  the  terms  are  the  gas  volume  fraction,  radial  gas 
velocity,  radial  liquid  velocity,  tangential  gas  velocity  and  tangential  liquid 
velocity  respectively.  The  data  format  is  free. 

3)  DELA,  AA 

To  determine  the  fraction  of  grid  size  in  which  the  injection  is  made. 

0 £ AA  _<  1.  This  is  used  to  define  the  location  of  jet.  The  region  of 
injection  could  cover  several  full  or  fractions  of  grid  sizes. 

4)  DERIVl,  Derivative 

To  get  the  first  derivative  of  a data  array  using  a center  difference  sch-  me 
except  the  two  end  points  in  which  three  points  near  the  boundary  are  used. 
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5)  DGCOEF,  Generalized  Coefficients 


To  calculate  the  added  mass,  drag  and  all  the  generalized  coefficients  (A  , 

di, 

A,  and  C. .).  This  is  the  heart  of  the  modeling, 
d iJ 

The  effective  coefficients  are  modeled  as: 


A = 
a 

"d  = 

* ''d2“d2 

^1  = 

= + 202/(1  + 3a^)) 

= a^02Pj_/(ct2  + 2o^/(l  + 302)) 

^dl 

2 

= 18 

^d2 

= 18  - 02/0.8)^82^) 

al 

na,/  na  na^  * 

= OL^  /(o^  + O2  ) 

^a2  = 

• 

= 1 - w 

al 

'^dl  = 

nd-,  nd  nd^ 

= O2  /(a^  + QI2  ) 

^^d2 

' 1 ■ “di 

6)  DWPDE, 

Partial  Differential  Equations. 

To  evaluate  the  values  of  the  increments  on  the  column  matrix  W from  the 
partial  differential  equations.  This  is  the  major  part  of  the  McCormack’s 
scheme.  In  each  complete  time  step  this  routine  will  have  to  be  called  twice 

7)  FNDDT,  At 

To  determine  the  suitable  time-step  size. 

8)  FSOFW,  Column  matrices  F and  S. 

To  determine  the  convective  matrix  F and  the  source  matrix  S. 
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9)  JET,  Injection. 


To  determine  the  momentum  source  due  to  the  jet  injection. 

10)  SIZES 

To  determine  the  gas  bubble  and  liquid  droplet  sizes.  In  the  model  the 
sizes  were  modeled  to  be  functions  only  of  the  volume  fraction,  i.e. 

d,  = d , a,  . 
k ok  k 

Different  models  for  size  distributions  could  be  easily  adopted  here. 

1 1 ) TAUOFW 

To  determine  the  stress  tensor  t and  its  derivative. 

12)  UOFW 

To  convert  the  column  matrix  W into  the  physical  independent  variables, 

such  as  a,  V , V^. 

r 0 

V.  Summary 

A computer  program  aimed  at  the  phase  separation  between  gas  and  liquid  at 
zero  gravity,  induced  by  vortex  motion,  is  developed.  The  vortex  motion  is 
created  by  fluid  injections.  The  computer  program  uses  a FORTRAN  77  based  code 
and  HP-1000  minicomputer.  It  is  flexible  and  accepts  various  input  parameters 
for  different  flow  conditions.  Other  interaction  effects  can  also  be  added  or 
modified  easily.  This  program  can  be  used  to  study  the  fluid  dynamical  behavior 
of  the  rotational  two-phase  fluids  in  a cylindrical  tank.  It  provides  a 
quick/easy  sensitivity  test  on  various  parameters  and  thus  provides  the  guidance 
for  the  design  and  use  of  actual  physical  systems  for  handling  two-phase  fluids. 
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Appendix  A 


Code  Listing 


iGLOM  T=00004  IS  ON  CR  T4 


USING  00136  BLKS  R=0000 


0001 

0002 

0D03 

0004 

0005 

0006 
0007 
0 0 08 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 
0 0 23 

0024 

0025 

0026 
0 0 27 
0028 
0 0 29 

0030 

0031 

0032 
0 0 33 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 
0 0 47 

0048 

0049 

0050 

0051 
0 0 52 

0053 

0054 
0 0 55 
0056 
0 0 57 
0058 


FTN77 

$EMA  /DATA/ , /UWUU/ , /COEFF/ , /SOURCE/ , /FANDS/ , /TAU/ 
$FILES  1 .2 

PROGRAM  GLUM< ,99), <860425.1537) 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS,  PROGRAM  WAS  DEVELOPED  TO  STUDY  THE  FLUID  DYNAMICAL  BEHAVIOR 
OF  A ROTATIONAL  TWO-PHASE  FLUIDS ( GAS/L lOUID ) IN  A CYLINDRICAL  TANK. 

THE  VORTEX  MOTIONS  ARE  ESTABLISHED  BY  TANGENTIAL  FLUID  INJECTION. 

THE  PROGRAM  WAS  DEVELOPED  ORIGINALLY  BY  T.T.  YEH  OF  NBS 
IT  WAS  BASED  ON  HP'S  FORTRAN  77(ANSI  77+MIL-STD-l 753 ) 

WHEN  WHO  WHAT 

8502XX  TTY  ZERO-G  FUEL  TRANSFER,  START-UP  STAGE. 

LAX-WENDROFF  2-STEP  SCHEME  (FULL  STEP  PREDICTION+CORRECTION) 
WITH  NUMERICAL  DAMPING  FACTOR ( NORMALLY  SET  TO  ZERO) 
CONSERVATION  FORM,  VAR  I ABLE < AUTO ) TIME  STEP 
REAL»8 

INTERFACIAL  FORCES:  DRAG, ADDED  MASS 

PUMP  CONDITION:  MOMENTUM  SOURCE  BUT  NO  MASS  SOURCE 
850715  TTY  GENERALIZED  EQUATIONS  AND  COEFF.  Cij 
851018  TTY  IN  ANSI  77  STANDARD<  WITH  A LITTLE  EXCEPTION  FOR 

TESTING  IN  HP-1000) 


C »»♦«*  INTERNAL  SUBROUTINES  ***** 

C DELA.,  DERIVl,  DGCOEF,  DWPDE  , FNDDT  , FSOF-W , 

C INIT,  JET,  SIZES,  TAUOFW  and  UOFW 

C MOST  OF  THE  LIST  OF  NOTATIONS  ARE  GIVEN  IN  SUBROUTINE .. INI T 

C 

CHARACTER  NAMR*16,  N0TES»72 
INTEGER  I,IOS,J,JTIME(5) , K , MM , NPRT , NT 
Z ,IW,LUP ,LUS,NG,NGM1 ,NGM2 

PARAMETER  <MM=101) 

REALi*8  BA(MM,2)  , DDT  , DT  , DTMAX  , DTMIN  , DTPRT  , P ZERO 
X ,DW1,  T,TMAX,TPRT, VJT,  VDR<2) 

Y ,RJ1 ,RJ2,TJ1 ,TJ2,QJ,VJ 

1 , ALMT,U, V, ALP ,P ,R , W , UP , WN , DW , RDP , RH , F,S 

4 ,BR,BRH,  RH0,MUEF,V18,NA,ND 

5 , DO, GAMMA,  DAMP, DR 

6 ,TRR,TRA,TAA,RTRR,RTRA,  C, CPA, CD 


COMMON 

Y /JETS/  RJl ,RJ2,TJ1 ,TJ2,QJ,VJ 

Z /CONTP/  IU,LUP,LUS,NG,NGM1 ,DAMP,DR 

1 /ALPLMT/  ALMT(2) 

2 /COEFF/  C<MM,2,2) ,CPA(MM,2) ,CD(MM,2) 

3 /DATA/  U(MM,2) ,V(MM,2) ,ALP(MM,2) ,P(MM) ,R(MM) 

4 /DRAGl/  RH0(4) ,MUEF(2) ,V18(2) ,NA,ND 

6 /FANDS/  F(MM,5) ,S<MM,5) 

7 /SOURCE/  BR(MM) ,BRH(MM) 

8 /TAU/  TRR(MM,2) ,TRA(MM,2) ,TAA<MM,2) ,RTRR(MM,2) ,RTRA(MM,2) 

9 /WWUW/  W<MM,5) ,UP(MM,5) ,UN(MM,5) ,DU(MM,5) ,RDP(MM) ,RH(MM> 

EQUIVALENCE  (WN,BA> 

C ******  RHO(l)  < RH0(2)  (i.e.  PHASE-1=GAS,  PHASE-2=L IQU I D ) 

C DTPRT  TIME  STEP  FOR  PRINTOUTCAND  PLOT) 

C 
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0 0 59 

0060 

0061 

0062 

0063 

0064 

0065 

0 0 66 

0067 

0068 

0 0 69 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0 0 77 

0078 

0079 

0080 

0081 

0082 

0083 

0084 

0085 

0086 

0087 

0 0 88 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0 0 96 

0097 

0098 

0099 

0100 

0101 

0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 

0111 

0112 

01  13 

0114 

01  15 

0116 

01  17 

0118 


LUP=1  ! LU  FOR  PRINTING  DEBUG  DATA(=1  TERMINAL) 

LUS=6  ! LU  FOR  STORING  DATA<*6  PRINTER) 

7 F0RMAT<2X,A,3I5) 

8 F0RMAT(2X,A,3(1PE12.4) ) 

URITE(1,7)  'Enter  lu  for  sawing  data.  D.F.=',LUS 
READd,*)  LUS 

C KEET  JOB  TIME  FOR  FUTURE  REFFERENCE 

C CALL  EXEC(11,JTIME,JTIME(D) 

IF (LUS  .NE.  1 .AND.  LUS  . NE . 6)  THEN 
C ********  Define  a file  nane  for  stroing  output  *»»****»*«*»*«»» 
URITEC 1 , ' (2A) ' ) 'Enter  FILE  NAME  for  sawing  data.' 

READ( 1 , ' (A) ' ) NAMR 
LUS=90 

OPEN(LUS,FILE=NAMR ,IOSTAT=IOS,STATUS='NEW' ,ERR=999) 

ENDIF 

URITE< 1 , ' (A) ' ) ' Enter  N0TES(<73  CHAR,)  for  the  job' 

READ( 1 , ' (A) ' ) NOTES 
URITE(LUS, ' (3H1  ,A)')  NOTES 

C URITE(LUS, ' (514) ' ) JTIME 

C To  set-up  the  initial  condition. 

CALL  INIT 
NGM2=NG-2 

VDR ( 1 >=2 . »MUEF ( 1 ) /DR*«2  ! For  deter«ining  tine  step 

UDR (2)*2 . »MUEF(2)/DR»»2 


T = 0 . ! 

TMAX=5, 

DTPRT=0 .2 

URITE(1,S)  'Enter  INITIAL  and 
READ(1,»)  T,TMAX 
URITE(1,8)  'Enter  TIME  STEP  f 
READ(1,«)  DTPRT 

DTMIN=1.0D-6  ! 

DTMAX=DTPRT 

NT=0  ! 

TPRT=T 

NPRT=0 

PZERO=0.D0  ! 

DT=DTMIN 


INITIAL  TIME 

FINAL  TIMES.  D . F . = ' , T , TMAX 
r output.  D.F.=',DTPRT 

SET  MINIMUM  TIME  STEP 

Ti«e  step  nuMber 

Pressure  at  center  core 


10  NT=NT+1 

CALL  UOFW(W,R,NG) 

PRINT  OUT  AT  SELECTED  TIME 
IF(T  .GE.  TPRT  .OR.  T . GT . TMAX)  THEN 
IF(NT  .GT.  1)  THEN 
NPRT=NPRT-M 

TPRT*TPRT+DTPRT»( 1 . +DNINT( ( T-TPRT ) /DTPRT ) ) 

FORMAT( IHl ,2( A5, 14, A5, 1PE9 . 3) ) 

URITE(LUP ,21 ) 'NP=' ,NPRT , 'T=' ,T, 'NT=' ,NT, '0T=' ,DT 
WRITE (LUS, 21 ) 'NP=' , NPRT , ' T= ' , T , 'NT=' ,NT, 'DT*' , DT 


WRITE (LUS, ' (A5,A6,5A1 1 ) ' ) ' J' , 'ALPl 


U1 


'U2' , 'VI 


V2' 
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0119 

0120 
0121 
0122 

0123 

0124 

0125 

0126 

0127 

0128 

0129 

0130 

0131 

0132 

0133 

0134 

0135 
01  36 

0137 

0138 

0139 

0140 

0141 

0142 

0143 

0144 

0145 

0146 

0147 

0148 

0149 

0150 

0151 

0152 

0153 

0154 

0155 

0156 

0157 

0158 

0159 

0160 
0161 
0162 

0163 

0164 

0165 

0166 

0167 

0168 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0178 


DO  30  J»1 ,NG 

30  URITE(LUS,  MI5,F6. 4,5(1  PEI  1 .3))') 

1 J , ALP (J,l),U(J,l),U<J,2),y(J,l),V(J,2),P<J) 

ENDIF 

ENDIF 

IF(  T .GT.  TMAX)  GOTO  9999 

C TO  SOLVE  THE  DIFFERENTIAL  EQUATIONS 
C USING  2 STEP  LAX-WENDROFF  SCHEME 

C MacCornack's  «0thod..  BACKWARD  PREDICTOR,  FORWARD  CORRECT 

C CENTER  DIFFERENCED  ON  TAU 

C FIRST  TO  GET  THE  SPECIAL  VARIABLES  AND  THEIR  SPACIAL  DERIVATIVES 

CALL  DGCOEF(ALP ,NG)  ! DRAG  AND  GENERIZED  COEFF . 

CALL  TAUOFW(MUEF,DR,R,NG)  ! STRESS 


VJT=VJ 

IF(  T .LT.  TJl  .OR.  T . GT . TJ2)  VJT=0 . ! NO  INJETCTION 

CALL  JET(BA,RHO,BR,V,VJT,NG) 

CALL  FSOFW(W,BA,R,NG) 

C DETERMINE  THE  TIME-STEP  SIZE 

IF(  NT  .GT.  1)  THEN 
CALL  FNDDT(DT,DR ,VDR,NG) 

DT=DT+5 . *DTMIN/1 0 . **NT 
IF(DT  .GE.  DTMIN)  THEN 
I=»DLOGl  0 (DT/DTMIN) 

DDT=DTMIN*1 0 . **I 
DT=DINT(DT/DDT+0 .001 )«DDT 
IF(  DT  .GT.  DTMAX)  DT=DTMAX 
ELSE 

WRITE(LUP, ' <5X,A,1PE13.3) ' ) 

1 'STOP  DUE  TO  TOO  SMALL  TIME  STEP.  DT-',DT 

GOTO  9999 
ENDIF 
ENDIF 

C BACKWARD  PREDICTOR 

CALL  DWPDE(DW,RDP,  DR , DT , RH , NGMl > ' INCREAMENT 

DO  40  J=1 ,NGM1 
DO  40  1=1 ,5 

40  WP( J, I >=0 . 5»(W< J+1 , I)+W( J, I ) )+DW( J, I ) ! BASE+INCREAM 

C PREDICTION  DATA  COMPLETED,  CONTINUE  FOR  CORRECTION 

CALL  UOFW(WP,RH,NGMl ) 

CALL  DGCOEF(ALP,NGMl ) 

CALL  TAUOFW(MUEF,DR,RH,NGMl ) 

CALL  JET{BA,RHO,BRH,V, VJT,NG) 

CALL  FSOFW(WP ,BA,RH,NGM1 ) 

C FORWARD  CORRECTION 

CALL  DWPDE(DW,P(2) , DR , DT , R ( 2 ) , NGM2 ) 

P(1 )=PZERO 
DO  50  J=1,NGM2 


! MOMENTUM  SOURCE 
! CONVECTIVE  -F  AND  SOURCE-S 

I INITIAL  INPLUSE  TREAMENT 
! ROUND  OFF  TIME  STEP 
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0179 

0180 

0181 

0182 

0183 

0184 

0185 

0186 

0187 

0188 

0189 

0190 

0191 

0192 

0193 

0194 

0195 

0196 

0197 

0198 

0199 

0200 

0201 

0202 

0203 

0204 

0205 

0206 

0207 

0208 

0209 

0210 

021 1 

0212 

0213 

0214 

0215 

0216 

0217 

0218 

0219 

0220 

0221 

0222 

0223 

0224 

0225 

0226 

0227 

0228 

0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 

0237 

0238 


50 


P < J+1 )=P< J)  + <0 .25»<RDP( J)+RDP ( J + 1 ) )+0 ,5*P( J+1 ) )/R ( J) 

DO  50  1=1,5 

UN<J+1,I)=.25»(WP<J+1,I)+UP<J,I))  + .5»<W<J+1 ,I)+DW<J,I)  ) 
P(NG)=P(NGM1 ) 

C 2ND  STEP  (PREDICTION+CORRECTION)  COMPLETED 

IF<  NT  .EQ,  1 ) GOTO  10 

C Estirtation  oP  initial  P coMpleted,  return  to  the  initial  condition 
C and  start  to  advance  the  program  in  time. 


»«»«««««««««««««««»«*« M »*«««««« 

C DATA  W AT  THE  NEW  TIME  STEP  COMPLETED 


C IMPOSED  B.C.  *6.4 

DW1=.5»(WN(2, 1 )+U<2, 1 ) )-0 . 

ALPd  ,1  ) = < ,5»<W<1  >5)+U(2,5))-DUl*DT/DR)/RH<l  ) 
IF(ALP(1,1)  .LT.  ALMT(D)  ALP  ( 1 , 1 ) =ALMT  ( 1 ) 
IF<ALP<1,1)  ,GT.  ALMT<2))  ALP ( 1 , 1 ) »ALMT ( 2 ) 
WN(1 ,5)=ALP( 1 , 1 )*R(1 ) 


DUl»0.-0 ,5»(UN<NGM1  .1)+U<NGM1  .1)) 


ALP<NG,1 )=( ,5»<U(NGM1 ,5)+U 
IF<ALP<NG,1)  .LT.  ALMT(1)> 
IF(ALP<NG,1)  .GT.  ALMT(2)) 
WN  < NG , 5 ) *ALP ( NG , 1 > »R ( NG ) 


NG,5) )-DUl»DT/DR )/RH(NGMl ) 
ALP(NG,1 )=ALMT(1 ) 

ALP<NG, 1 )=ALMT(2) 


UN(1,1)=0.  ■ ! NO  RADIAL  VEL . 

WN(1 ,2)=0 . 

UN(NG, 1 )=0 . 

UN<NG,2)=0 . 

UN(1,3)=UN<2,3)*WN(1 ,5)/UN(2,5)»R<l)/R(2)  ! OMEGA=CON. 

WN  < 1 , 4 ) =WN  < 2 , 4 ) * ( R ( 1 ) -WN  <1,5))/(R(2) -WN (2,5))*R(1)/R(2) 
IF<IW  .EQ.  0)  THEN 

WN(NG,3)=WN(NGM1 , 3 ) *WN< NG , 5 ) /WN ( NGMl , 5 ) *R < NG ) /R < NGMl ) 
WN(NG,4)=WN(NGM1 ,4)*<R<NG)-WN(NG,5) )/(R(NGMl )-WN(NGMl ,5) ) 
1 *R<NG)/R<NGM1 ) 


ELSE 

WN(NG,3)»0.  ! NON-SLIP  AT  WALL 

WN(NG,4)=0 . 

ENDIF 


C ARTIFICIAL  DAMPING 
DO  60  1*1,5 

W(1  ,I)  = <1  ,-DAMP)*WN(l , I )+DAMP*WN(2, 1 ) 
W(NG,I)=(1 .-DAMP)»WN<NG,I)+DAMP»WN(NGM1 ,1) 
DO  60  J*2,NGM1 


60 

W(J,I)  = <1  , 

-DAMP ) *WN<  J , I ) +DAMP*  < WN( J-1 

,I)+WN(J+1,I)-WN(J,I)) 

C 

SOLUTION  FOR 

THIS  TIME  STEP  COMPLETED 

90 

T-T+DT 
GOTO  10 

! UPDATE 

TIME 

AND  CONTINUE  TO  THE  NEXT 

999 

WRITE(LUP,7) 

WRITE(LUP,7) 

WRITE<LUP,7) 

'OPEN  FILE  FAILED  ON 
NAMR 

' lOSTAT*' , lOS 

FILE: 

/ 

9999 

CONTINUE 
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0239 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

0247 

0248 

0249 

0250 

0251 

0232 

0253 

0254 

0255 

0256 

0257 

0250 

0259 

0260 

0261 

0262 

0263 

0264 

0265 

0266 

0267 

0268 

0269 

0270 

0271 

0272 

0273 

0274 

0275 

0276 

0277 

0278 

0279 

0280 

0281 

0282 

0283 

0284 

0285 

0286 

0287 

0288 

0289 

0290 

0291 

0292 

0293 

0294 

0295 

0296 

0297 

0298 


C CALL  EXECCll ,JTIME,JTIME<1 ) ) 

C URITE(LUS, ' (35X,5I4) ' ) JTIME 

CLOSE (LUS) 

END 


»«»««*««««**»«»*«»«»» »»«••«» 

$EMA  /DATA/, /WUWU/, /SOURCE/ 

SUBROUTINE  INIT  ,< 860425 . 1 537 > 

C TO  SET-UP  THE  INITIAL  CONDITIONS 
INTEGER  I,IOS,ITLOC,IVTX, 

Z ,ITIME,IW,LUP,LUS,NG,NGM1 

PARAMETER  (MM*101) 

CHARACTER  NAMR*16 
REAL»8  ALMT,  W , WP , WN , DU , RDP , RH 
1 ,U,V,ALP,P,R,  F,S 

4 ,RHO,MUEF,V18,NA,ND,  DO, GAMMA 

7 ,BR,BRH,  DAMP, DR 

Y ,RJ1 ,RJ2,TJ1 ,TJ2,QJ,VJ 

X ,A2,DELA,DEN1D2,D1 ,D2,PI ,RE,RMIN,RTANK 

X ,DS,L3,VS,TS,PS 

X ,ALP10, OMEGA, RPEAK,RJB,VPEAK 

X ,EUF<2) ,VT,MU<2) ,MU1D2 

COMMON 

1 /ALPLMT/  ALMT(2) 

Y /JETS/  RJl  ,RJ2-,TJ1  ,TJ2,QJ,UJ 

Z /CONTP/  IW,LUP,LUS,NG,NGM1 ,DAMP,DR 

3 /DATA/  U(MM,2) ,V(MM,2) ,ALP<MM,2) ,P(MM) ,R(MM) 

4 /DRAGl/  RH0<4) ,MUEF(2> ,yi8(2) ,NA,ND 

3 /DSIZE/  D0(2) ,GAMMA(2) 

7 ./SOURCE/  BR  (MM)  ,BRH(MM) 

9 /UUWU/  U(MM,5) ,UP(MM,5) ,UN(MM,5) ,DW(MM,5) ,RDP(MM) ,RH(MM) 

DATA  PI/3. 141596D0/ 


C »*«***  NOTES: 

C RHO 

C ALMT 

C DAMP 

C DS,LS,yS,TS 

C RE 

C RTANK 

C IVTX 

C 

C RPEAK,VPEAK 

C OMEGA 

C EVF,MUEF 

C DO, GAMMA 

C RJ1,RJ2,TJ1, 

C TJ2,QJ,VJ 

C 


PHASE-1=»GAS,  PHASE-2»LIQUID 
DENSITY,  RHO(l)  < RH0<2> 

LIMIT  VALUES  FOR  ALPl,  ALM ( 1 ) < ALP  1 < ALM ( 2 ) 
NUMERICAL  DAMPING  FACTOR,  NORMALLY  =0. 
DENSITY, LENGTH, VELOCITY  AND  TIME  SCALES 
REYNOLDS  *=Vs»R  TANK  *»RHO  ( 2 ) /MU  < 2 ) 

TANK  RADIUS 

TYPE  OF  INITIAL  FLOW. 

( 0=Rest , 1 *P ur«  r o Ta T i on , 2»H . 0 , 3“GIT  Vortex) 

VORTEX  PARAMETERS 

PURE  ROTATION.  V=»OMEGA*r 

EFFECTIVE  VISCOSI TY  , MUEF=* ( 1 t-EVP  ) »MU 

DIA.  PARAMETERS:  D=DO*ALP »»GAMMA 

TO  DEFINE  JET  SIZE,  PUMPING  TIME, 

VOLUME  FLOW  RATE  AND  INJECTION  MEAN  SPEED 


5 F0RMAT(2X,A,2(X,A) ) 

7 F0RMAT(2X,A,3I5) 

8 F0RMAT<2X,A,3(1PE12.4) ) 

9 F0RMAT(A25,2(1PE15.4) ) 

92  FORMAT(X,7( IPEl 1 . 4) ) 
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0299 

0300  C DEFINE  THE  PARAMETERS  FOR  THE  PROBLEM. 

0301 

0302  RTANK=1.D0  ! COULD  BE  SET  TO  1 (M) 

0303  RH0(2)=1 . OOOD+3  ! " (KG/M»»3) 

0304  MU(2)=1 .514D-3  ! “ (KG/M-S) 

0305 

0306  C WRITE(1,8)  'ENTER  RTANK<M)  OR  DEFAULT  ',RTANK 

0307  C READ(1,»)  RTANK 

0308 

0309  C The  values  of  RTANK , RHO ( 2 ) , and  MU<2)  could  all  be  set  to  1,  since 

0310  C the  length  and  density  scales  are  based  on  RTANK, and  RH0<2)  and  the 

0311  C value  of  the  viscosity  MU(2)  can  be  coHbined  into  and  specified  by  the 

0312  C Reynolds  nunber  RE.  Thus  all  character istic  sea les ( LS , VS , TS,and  DS  ) 

0313  C are  fixed  after  RE  is  given, 

0314 

0315  RE=1.0D5  ! LS««2/ < NU»TS ) =LS*VS/NU 

0316  WRITE(1,8)  'enter  Reynolds  no., RE.  D.F.=',RE 

0317  READ(1,»)  RE 

0318  LS=RTANK  ! LENGTH  SCALE  (M) 

0319  DS=RH0<2)  ! DENSITY  SCALE  (KG/M»»3) 

0320  VS=RE»MU(2)/RH0(2)/LS  ! VELOCITY  SCALECM/S) 

0321  TS=LS/VS  ! TIME  SCALE  (S) 

0322  PS=DS»VS**2  ! PRESSURE  SCALE 

0323 

0324 

0325  C AFTER  THIS  POINT  ALL  VARIABLES  ARE  BAESED  ON  THE  CHARA.  SCALES 

0326  C i.e.  ALL  VARIABLES  ARE  DIMENSIONLESS 

0327 

0328 


0329 

DEN1D2=1 .293D0/1 .000D3 

! 

D1/D2 

0330 

0331 

0332 

MU1D2=1 .71D-5/1 .514D-3 
RH0(2)=RH0(2)/DS 

1 

MU1/MU2 

0333 

0334 

0335 

0336 

MU<2)=MU(2)/(DS»LS»VS) 
RHO( 1 )=DEN1D2*RH0(2) 
MU< 1 )=MU1D2»MU(2) 

! 

= 1/RE 

0337 

ALMT(1)=0. 00010 

1 

MIN.  OF  ALPl 

0338 

ALMT(2)=0 . 99999 

1 

MAX.  OF  ALPl 

0339 

DO( 1 )=1 .D-2/LS 

! 

GAS  DIAMETER  at  ALP1=»1 

0340 

D0(2)=l .D-2/LS 

1 

LIQUID  DIAMETER  at  ALP2=1 

0341 

0342 

GAMMAd  )=2.D-1 
GAMMA<2)=2. D-1 

j 

0343 

EVF(2)=1 . D3 

1 

TURB.+PHASE-DISPERSION  EFFECTS 

0344 

EVF( 1 )=DEN1D2/MU1D2»EVF(2 

) ! 

MODEL 

0 345 

DAMP=0 . 0 

I 

NUMERICAL  DAMPING  FACTOR < e . g .=. 2 ) 

0346 

NA=4. 

! 

WEIGHTING  EXP.  FOR  ADM 

0347 

0348 

ND=4. 

! 

WEIGHTING  EXP.  FOR  DRAG 

0349 

WRITE(1,8)  'Enter  DENSITY 

and 

VISCOSITY  ratios. ' 

0350 

0351 

WRITE(1,8)  'D.F.*' ,DEN1D2 
READ<1,»)  DEN1D2,MU1D2 

,MU1D2 

0352 

0353  URITE(1,8)  'Enter  BASE  DIAMETERS:  D01,D02' 

0354  WRITE(1,8)  'D.F.=',DO 

0355  READ<1,»)  DO 

0356 

0357  URITE(1,8)  'Enter  SIZE  EXPONENT:  GAMMAl , GAMMA2' 

0358  WRITE(1,8)  'D.  F.  =»', GAMMA 
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0359 

0360 

0361 

036x2 

0363 

0364 

0365 

0366 

0367 

0368 

0369 

0370 

0371 

0372 

0373 

0374 

0375 

0376 

0377 

0378 

0379 

0380 

0381 

0382 

0383 

0384 

0385 

0386 

0387 

0388 

0389 

0390 

0391 

0392 

0393 

0394 

0395 

0396 

0397 

0398 

0399 

0400 

0401 

0402 

0403 

0404 

0405 

0406 

0407 

0408 

0409 

0410 

041 1 

0412 

0413 

0414 

0415 

0416 

0417 

0418 


READ<1,»)  GAMMA 

WRITE(1,8)  'Enter  weighting  exponent:  NA,ND.  D.F,=',NA,ND 
READU,*)  NA,ND 

WRITE(1,8)  'Enter  GAS  VOLUME  FRACTION  1 ini ts : ALMT 1 , ALMT2 . ' 
WRITE(1,8)  'D.F,*',ALMT 
READ<1,»)  ALMT 

URITE(1,8)  'Enter  eddy  wisicosity  factor.  D.F.=',EVF 
READd,*)  EVF 


IW=1 

WRITE(1,7)  'Enter  wall  condition , l=nonslip , 0=slip . D.F.',IW 
READd, »)  IW 

WRITE(1,8)  'Enter  nunerical  damping  factor.  D.F.=',DAMP 
READd,*)  DAMP 

DO  10  K=l,2 
V18(K)=18.*MU(K) 

10  MUEF(K)=MU<K)»d .+EVF(K) ) ! EFFECTIVE  VISCOSITY  FOR  STRESS 

RH0(3)=RH0d  )»RH0<2) 

RH0(4)=RH0d  )-RH0(2) 

RMIN=0.1  ! MINIMUM  FLOW  RADIUS  IN  THE  TANK 

NG=101  ! ♦ OF  GRID  POINTS  USED 

NGM1=NG-1 
DR  = d . -RMIN)/NGM1 

C Initial  c leann ing-up . 

DO  15  J=1,MM 
DO  15  K=1 ,6 
W( J,K )=0 . DO 
15  U<J,K)=0.D0 

C MOMENTUM  SOURCE,  JET  CONDITIONS 

RJ1=8.5D-1 
RJ2*9.5D-1 

VJ=10.  ! TANGENTIAL  INJECTION  SPEED 

TJ1=0 . 

TJ2=«10  . 

WRITEd,8)  'Enter  JET  SIZE  defined  by  RJ1,RJ2.  D . F . - ' , R J1  , R J2 
READd,*)  RJ1,RJ2 

WRITE(1,8)  'Enter  INJECTION  SPEED  AND  TIME  RANGE,  VJ,T1,T2' 
URITEd,8)  'D.F.=' ,VJ,TJ1  ,TJ2 
READd,*)  VJ,TJ1,TJ2 

QJ=<RJ2-RJ1 )*VJ  ! JET  VOLUME  FLOU  RATE 

DO  20  J=t ,NG 
R( J)=RMIN+< J-1 )*DR 
RH<  J)=R(J)t-0 .5*DR 

BR < J)=DELA( (R < J) ) ,DR ,R J1 ,RJ2)/(2 . *PI ) ' JET  DISTRIBUTION 

20  BRH(J)=DELA( (RH(J) ) ,DR,RJ1 ,RJ2)/(2. »PI)  I PER  RADIAN 


C SETUP  INITIAL  CONDITIONS 

IVTX=0 
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0419 

0420 

0421 

0422 

0423 

0424 

0425 

0426 

0427 

0428 

0429 

0430 

0431 

0432 

0433 

0434 

0435 

0436 

0437 

0438 

0439 

0440 

0441 

0442 

0443 

0444 

0445 

0446 

0447 

0448 

0449 

0450 

0451 

0452 

0453 

0454 

0455 

0456 

0457 

0458 

0459 

0460 

0461 

0462 

0463 

0464 

0465 

0466 

0467 

0468 

0469 

0470 

0471 

0472 

0473 

0474 

0475 

0476 

0477 

0478 


OMEGA=0 . 
yPEAK=0 . 
RPEAKsRMIN 


NAMR='SiHple  uortex' 

URITE(1,5)  'Enter  data  FILE  NAME  for  initial  cond.,  if  any;' 
WRITE(1,5)  'D.F,=',NAMR 
READ(1,'(A)')  NAMR 

IF<  NAMR  .NE.  ','  .AND.  NAMR  . NE . 'Simple  vortex')  THEN 
C INITIAL  CONDITION  FROM  A GIVEN  FILE  NAMR. 

OPEN ( 99 , F ILE=NAMR , IOSTAT=IOS , STATUS= ' OLD ' , ERR=299 ) 

DO  25  J»1,NG  ! INITIAL  VALUES  FROM  FILE  NAMR 

25  READ<99,*)  K , ALP < J , 1 ) , U < J , 1 ) , U ( J , 2 ) , V ( J , 1 ) , V C J , 2 ) 

CLOSE (99) 

ELSE 

C TO  DEFINE  INITIAL  CONDITION. 

ALP10=2.5D-1  ! INITIAL  GAS  VOL.  FRACTION 

WRITE(1,8)  'Enter  initial  value  of  alpl.  D.F.='.ALP10 
READ(1,»)  ALPIO 

URITE(1,7)  'Enter  type  of  vortex:  0=At  rest,l=pure  rotation' 
URITE(1,7)  '2=H.O. ,3=GIT.  D.F.=',IVTX 
READ(1,»)  IVTX 
IFdVTX  .GT.  0)  THEN 
IF (IVTX  .GT.  1)  THEN 

URITE(1,3)  'Enter  PEAK  VEL . and  LOCATION  for  classic  vortex' 
WRITE(1,8)  'D.F.=' ,VPEAK,RPEAK 
READ(1,»)  VPEAK,RPEAK 

IF(RPEAK  .LE.  0.)  RPEAK=1 . ! SINGULAR  AT  ZERO 

ELSE 

WRITE(1,8)  'Enter  CIRCULAR  SPEED(rad . /unit  time).D.F.=' 

1 , OMEGA 

READ(1,»)  OMEGA 
ENDIF 
ENDIF 

DO  30  J=1 ,NG 
ALP( J, 1 )=ALP1 0 

VT=OMEGA»R( J)  ! PURE  ROTATION 

IFdVTX  .GT.  0)  THEN 
RJB=R( J)/RPEAK 

IF(  IVTX  .EQ.  1)  THEN  ! H.O.  VORTEX 

VT=VT+1 . 398»VPEAK/RJB*( 1 ,-DEXP(-l . 25643»R JB*»2 ) ) 

ELSE  I G.I.T.  VORTEX 

VT=VT+VPEAK»RJB»DEXP ( ( 1 . -RJB«*2)/2. ) 

ENDIF 

ENDIF 

DO  30  K=*l  ,2 
IFdU  .EO.  0)  THEN 

V(J,K)=»VT  ! NO  WALL 

ELSE 

V( J,K )=VT»( 1 . -R d) )»»0 . 1 ! BOUNDARY  LAYER 

ENDIF 

30  U(J,K)=0.  ! NO  RADIAL  VEL. 

ENDIF 

DO  40  J*1,NG  ! FORM  U FOR  NUMERICAL  CAL. 
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0479 

0480 

0481 

0482 

0483 

0484 

0485 

0486 

0487 

0488 

0489 

0490 

0491 

0492 

0493 

0494 

0495 

0496 

0497 

0498 

0499 

0500 

0501 

0502 

0503 

0504 

0505 

0506 

0507 

0508 

0509 
051  0 

051 1 

0512 

0513 

0514 

0515 

0516 

0517 

0518 

0519 

0520 

0521 

0522 

0523 

0524 

0525 

0526 

0527 

0528 

0529 

0530 

0531 

0532 

0533 

0534 

0535 

0536 

0537 

0538 


ALP<J,2)=1 .-ALP(J,1 > 

P(J)=O.DO 

W(J,5)=R(J)»ALP(J, 1 ) 

DO  40  K = 1 ,2 

W ( J , K ) =ALP  < J , K ) *U  < J , K ) »R  < J ) 

40  W(J,K+2)=ALP(J,K)»V(J,K)»R(J) 


PRINTOUT  PARAMETERS 


URITE 

WRITE 

WRITE 

WRITE 

WRITE 

WRITE 

WRITE 

WRITE 

URITE 

WRITE 

WRITE 

URITE 


(LUS,5) 
<LUS,5) 
<LUS,9) 
(LUS,9) 
<LUS,9) 
(LUS,9) 
(LUS,9) 
(LUS,*) 
<LUS,9) 
(LUS, 9) 
(LUS, 9) 
(LUS, 9) 


'INITIAL  CONDITION  FILE:',NAMR 
'••DIMENSION  UNITS  ARE  IN  MKS^»' 
'DENSITY  SCALE(kg/M3) ' ,DS 
'LENGTH  SCALE=RTANK, (m)' ,LS 
'VELOCITY  SCALE(m/5) ' ,VS 
'TIME  SACLE(s) ' ,TS 
'PRESSURE  SCALE(Pa)' ,PS 

'Reynolds  nuMber,  Re', RE 
'Jet  size,  RJl ,RJ2' ,RJ1 ,RJ2 
'Tangential  jet,  QJ,VJ',QJ,VJ 
'Injection  t i«e , T J1 , TJ2 ' , TJl ,TJ2 


URITE(LUS, ' (/33X, "PHASE-l 8X , “PHASE-2 ") ' ) 

WRITE (LUS, 9)  'Density' ,RHO( 1 ) ,RH0(2) 

URITE(LUS,9)  'Viscosity ', MU 
WR ITE( LUS , 9 ) 'Eddy  viscosity  factor',EVF 
URITE(LUS,9)  'Base  dia . ' ,DO 
WRITE(LUS,9)  'Size  exp.', GAMMA 

WRITE(LUS,9)  'Phase  1 ini t s ' , ALMT ( 1 ) , 1 . -ALMT ( 2 ) 

WRITE(LUS,^) 

WRITE(LUS,*) 

WRITE(LUS, ' ( “ OTHER  CONSTANS : lU , I VTX , NA , ND , DAMP , VPEAK , R PEAK ' 
1 ■*,0MEGA,D1/D2,MU1/MU2")  ' ) 

WRITE(LUS, ' (2I3,2F4. 0 ,F5. 2) ' ) lU , I VTX , NA , ND , DAMP 
WRITE(LUS,92)  VPEAK , RPEAK , OMEGA , DENI D2 , MUl D2 
URITE(LUS,») 

WRITE(LUS, ' (II 0,F10 .2) ' ) NG,RMIN 


RETURN 

299  WRITE(LUP,5)  'OPEN  FILE  FAILED  ON  INPUT  FILEi',NAMR 
WRITE(LUP,7>  'IOSTAT=' , lOS 
STOP  111 
END 


••»»»»*»**»**»»••••••»•••••*•••*»«•••*••*•*••*•••*•*****»***»***• 

REALMS  FUNCTION  DELA ( R , DR , RJ 1 , RJ2 ),( 860425 . 1 537 > 

C TO  DETERMINE  THE  EFFECTIVE  NOZZLE  SIZE  AT  EACH  GRID  LOCATION 
C THE  SIZE  IS  IN  THE  FRACTION  OF  GRID  SIZE  DR  (i.e.  0<DELA<1) 

REAL*8  R ,DR,RJ1 ,RJ2,  R1,R2 

R1=R-0.5^DR 
R2=R1+DR 
DELA=0 .DO 

IF(R1  .GE.  RJ2  .OR.  R2  . LE . RJl)  RETURN 

IF(R1  ,LT.  RJl)  R1=RJ1 

IF(R2  ,GT.  RJ2)  R2=RJ2 

DELA=(R2-R1 )/DR 

RETURN 
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0539 

0540 

0541 

0542 

0543 

0544 

0545 

0546 

0547 

0548 

0549 

0550 

0551 

0552 

0553 

0554 

0555 

0556 

0557 

0558 

0559 

0560 

0561 

0562 

0563 

0564 

0565 

0566 

0567 

0568 

0569 

0570 

0571 

0572 

0573 

0574 

0575 

0576 

0577 

0578 

0579 

0580 

0581 

0582 

0583 

0584 

0585 

0586 

0587 

0588 

0589 

0590 

0591 

0592 

0593 

0594 

0595 

0596 

0597 

0598 


END 


SUBROUTINE  DERI VI (Y,DY,DX,N2) , <860425. 1537> 

C GET  1ST  DERIVATIVE,  USING  CENTERED  DIFFERENCE 

REAL»8  DX,Y(1 ),DY(1) ,C 
EMA  Y,DY 

C=5 .D-l/DX 
DO  10  J=2,N2-1 

10  DY< J)=C»<Y( J+1 )-Y< J-1 ) ) 

*V  CALL  DWSUB(Y(3) ,1 , Y, 1 ,DY<2) , 1 ,N2-2) 

DY< 1 )=<Y<2)-Y( 1 ) )/DX  ! BASED  ON  3-END  PTS 

DY(N2)=(Y(N2)-Y(N2-1 > )/DX 
DY(1 )=2.*DY(1 )-DY<2) 

DY(N2)=2. »DY(N2)-DY(N2-1  ) 

*V  CALL  DUSMY<5. D-l/DX, DY, 1 ,DY, 1 ,N2) 

RETURN 

END 


♦EMA  /COEFF/  . 

SUBROUTINE  DGCOEF < ALP , N2 ),< 860425 . 1 537 > 

C CALCULATE  THE  DRAG,  ADDED  MASS  AND  GENERIZED  COEFF, 

INTEGER  J,MM,N2 
PARAMETER  (MM=«101) 

REAL»8  ALP(MM,2) 

EMA  ALP 

REAL»8  C, CPA, CD,  RHO , MUEF , V 1 8 , NA , ND 
X ,AA,AA1 ,AA2,AD,AD1 ,AD2,A1 ,A2,A12,DA12,DB2,  D 1 , D2 , WTl , UT2 , X 

COMMON 

1 /COEFF/  C(MM,2,2) ,CPA(MM,2) ,CD(MM,2) 

4 /DRAGl/  RH0(4) ,MUEF(2) ,V18(2) ,NA,ND 

DO  50  J=1 ,N2 
A1=ALP(J, 1 ) 

A2=ALP( J,2) 

A12=A1»A2 

C TO  GET  DRAG  COEFF.  AD 

CALL  SIZESCDl ,D2,A1  ) 

AD=V18( 2)»Al/( A2*D1*D1 ) ! =AD1  IF  A2>.78 

IF(  A2  .LT.  ,78)  THEN 
AD2=»V18(1  )»A2/<  (1  ,-A2/.8)*D2)**2 
IF(AD2  .LT,  AD)  THEN 
WT1=A2*»ND 
WT2=A1»»ND 

AD=<AD*UT1+AD2*UT2)/(WT1+WT2) 

ENDIF 

ENDIF 

C ADDED  MASS  COEFF.  AA 
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0599 

0600 

0601 

0602 

0603 

0604 

0605 

0606 

0607 

0608 

0609 

0610 

061 1 

0612 

0613 

0614 

0615 

0616 

0617 

0618 

0619 

0620 

0621 

0622 

0623 

0624 

0625 

0626 

0627 

0628 

0629 

0630 

0631 

0632 

0633 

0634 

0635 

0636 

0637 

0638 

0639 

0640 

0641 

0642 

0643 

0644 

0645 

0646 

0647 

0648 

0649 

0650 

0651 

0652 

0653 

0654 

0655 

0656 

0657 

0658 


AAl=A12»RH0(2)/(Al+A2/( .5+1 ,5*A1 ) ) 
AA2=A12*RH0(1 )/(Al/< .5+1 .5*A2)+A2) 

WT1=A2»*NA 

WT2=A1»»NA 

AA= ( AA 1 «WT1 +AA2»UT2 ) / ( UT 1 +UT2 ) 

C THE  GENERIZED  COEFF . CPA, C, AND  CD 

DB2=A12»RH0(3)+AA»(RH0( 1 ) »A 1 +RHO < 2 ) »A2 ) 
CPA< J, 1 )=A1«<A12*RH0(2)+AA)/DB2 
CPA< J,2)=A2»(A12*RH0(1 )+AA)/DB2 
C( J,1 ,2)*AA/DB2 
C< J,2, 1 )=C<J,1 ,2) 

C(J,1 ,1 )=A2*RH0(2)/DB2+C<J,1 ,2) 

C( J,2,2)=A1»RHQ(1 )/DB2+C( J,2, 1 ) 

CD< J, 1 )=- A2*RHQ(2)»AD/DB2 
CD< J,2)=A1»RH0< 1)»AD/DB2 

50  CONTINUE 

RETURN 

END 


$EMA  /DATA/ , /FANDS/ , /TAU/ , /COEFF/ 

SUBROUTINE  DUPDE ( DU , RDP , DR , DT , RR , N2 ),< 860425 . 1 537 > 

C 

C TO  GET  DU  OF  THE  PDEs 
C 

INTEGER  J, JPl ,K,KP2,MM,N2 
PARAMETER  (MM=»101) 

REAL»8  DW<MM,5) ,RDP(MM) , DR,DT,RR(MM) 

EMA  DW,RDP,RR 

C NOTES  1 COEFF,  C =»C*ALP  WHEN  THIS  IS  CALLED 

REAL*8  C, CPA, CD,  U,V,ALP,P,R,  F,S 

8 ,TRR,TRA,TAA,RTRR,RTRA 

4 ,RH0,MUEF,U18,NA,ND 

X ,ALP1 ,ALP2,CP1 ,CP2,DTDR,DU1 ,DU1DT,G1 , G2 , HDT , UT , W J 1 ,UJ3,UJ4 

COMMON 

1 /COEFF/  C(MM,2,2) ,CPA(MM,2) ,CD(MM,2) 

3 /DATA/  U(MM,2) ,U(MM,2) ,ALP <MM,2) ,P(MM) ,R (MM) 

4 /DRAGl/  RH0(4) ,MUEF(2) , V18<2) ,NA,ND 
6 /FANDS/  F(MM,5) ,S(MM,5) 

3 /TAU/  TRR(MM,2) ,TRA(MM,2) ,TAA(MM,2) ,RTRR(MM,2) ,RTRA(MM,2) 

DTDR=DT/DR 
HDT=0 . 5*DT 


DO  10  J=l,N2+l  ! CHANGE  C TO  ALP*C 

C( J, 1 , 1 )=ALP ( J, 1 )«C( J, 1 , 1 ) 

C< J,2, 1 )=ALP( J,2)»C( J,2, 1 ) 

C( J, 1 ,2)=ALP ( J, 1 )»C< J, 1 ,2) 

10  C(J,2,2)=ALP(J,2)*C< J,2,2) 


DO  20  J»1 ,N2 
JP1=J+1 

DU<  J,5)=DTDR)*(-F(  JPl  ,5)+F(J,5)  )+HDT»(S(  JPl  ,5)+S(  J,5)  > 
DO  25  K=1 ,2 
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0659 

0660 

0661 

0662 

0663 

0664 

0665 

0666 

0667 

0668 

0669 

0670 

0671 

0672 

0673 

0674 

0675 

0676 

0677 

0678 

0679 

0680 

0681 

0682 

0683 

0684 

0685 

0686 

0687 

0688 

0689 

0690 

0691 

0692 

0693 

0694 

0695 

0696 

0697 

0698 

0699 

0700 

0701 

0702 

0703 

0704 

0705 

0706 

0707 

0708 

0709 

071  0 

0711 

0712 

0713 

0714 

0715 

0716 

0717 

0718 


G1=.5*(C(JP1 ,K,1)+C(J,K,1)  ) 

G2=  . 5*  ( C ( JP 1 , K , 2 ) +C  ( J , K , 2 ) ) 

DU( J,K)=DTDR»(-F< JPl , K ) +F ( J , K ) +G1 » < RTRR < j P 1 , 1) -R TR R < J , 1 ) ) 

1 +G2»(RTRR(JP1 ,2)-RTRR( J,2) ) ) 

2 +HDT»(S< JPl ,K)+S( J,K) ) 

KP2=K+2 

25  DW< J,KP2)=DTDR*(-F( JPl ,KP2)+F(J,KP2) 

1 +G1»(RTRA< JPl , 1 )-RTRA( J,  1) ) 

2 +G2*(RTRA( JPl ,2)-RTRA( J,2) ) ) 

3 +HDT*(S( JPl ,KP2)+S( J,KP2) ) 

C DP  FOR  PRESSURE  CORRECTION 

CP1=0 .5*(CPA( J, 1 )+CPA( JPl  , 1 ) ) 

CP2=0 .5*(CPA< J,2)+CPA< JPl ,2) ) 

IF<  -DU<J,1)  .GT.  DW(J,2))  DW < J , 1 ) =-DU < J , 2 ) ! DP>=0 

RDP ( J)=(DU( J, 1 )+DU( J,2) )/(CPl+CP2) 

DU( J,1 )=DW< J,1 )-CPl»RDP( J) 

DU( J,2)=-DW< J,1 ) 

20  RDP ( J)=RDP ( J)/DT 

RETURN 

END 


*EMA  /COEFF/,/DATA/,/UWWU/ 

SUBROUTINE  FNDDT ( DT , DR , UDR , NG ),< 860425 , 1 537 > 

C DETERMINE  THE  TIME-STEP  SIZE 


INTEGER  I , J,LUP,MM,NG 
PARAMETER  (MM=101) 

REAL»8  DT,DR,VDR(2) 

REAL»8  C, CPA, CD,  RHO , MUEF , V 1 8 , NA , ND 
9 ,y,WP,UN,DW,RDP,RH,  U,V,ALP,P,R 

X ,DUM1,DUM2 


COMMON 

2 /COEFF/  C(MM,2 

3 /DATA/  U(MM,2) 

4 /DRAGl/  RH0<4) 
9 /UUUU/  U<MM,5 

DATA  LUP/1/ 


,2) ,CPA(MM,2) ,CD(MM,2) 

,U(MM,2) , ALP (MM, 2) ,P(MM) ,R(MM) 

,MUEF(2) ,V18<2) ,NA,ND 

) ,WP(MM,5) ,WN(MM,5) ,DU(MM,5) ,RDP(MM) ,RH(MM) 


DUM1=0 . 

DO  10  J=1,NG 
DUM2=-CD( J, 1 ) 

1 t-DAFS(U(  J,  1 ) )/DR 

2 +UDR( 1 )«ALP( J, 1 )*C( J, 1 , 1 ) 

3 +UDR(2)*ALP( J,2)*C( J,  1 ,2) 
IFCDUMl  .LT.  DUM2)  DUM1=DUM2 
DUM2=CD( J,2)+DABS(U( J,2) )/DR 

1 +UDR<2)*ALP(J,2)»C(J,2,2) 

2 >UDR(1 )»ALP(J,1)»C(J,2,1  ) 
IF(DUM1  .LT,  DUM2)  DUM1=DUM2 

10  CONTINUE 


' 1/St 

! CONVECTIVE 
' VDR=2»MUEF/DR**2 
! CROSS  VISICOSITY 


C FIND  THE  MAXIMUN  OF  DU 
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0719 

0720 

0721 

0722 

0723 

0724 

0725 

0726 

0727 

0728 

0729 

0730 

0731 

0732 

0733 

0734 

0735 

0736 

0737 

0738 

0739 

0740 

0741 

0742 

0743 

0744 

0745 

0746 

0747 

0748 

0749 

0750 

0751 

0752 

0753 

0754 

0755 

0756 

0757 

0758 

0759 

0760 

0761 

0762 

0763 

0764 

0765 

0766 

0767 

0768 

0769 

0770 

0771 

0772 

0773 

0774 

0775 

0776 

0777 

0778 


CALL  DUMAXC I ,DU, 1 ,NG)  ! VECTOR  OPERATION 

»V  DUM1=DABS(DW(I, 1 ) ) 

»V  CALL  DUMAX<J,DW(1 ,2) ,1 ,NG) 

*V  DUM2=DABS<DW< J,2) ) 

DT=1 .DO/DUMl 

RETURN 

END 

«»««««««« *««»««««*»«««««»»«*«*«»«««»«««»«««««««««««»**««««« «««»»««»* 
$EMA  /DATA/ , /FANDS/ , /TAU/ , /COEFF/ 

SUBROUTINE  FSOFU ( U , BA , RR , N2 ) ,<860425. 1537> 

C 

C CALCULATE  THE  CONVECTIVE-F  AND  SOURCE-S  TERMS 

C 

INTEGER  J,K ,KP2,L,MM,N2 
PARAMETER  (MM=101) 

REAL*8  U,U,ALP,P,R,  F,S,  RHO , MUEF , V 1 8 , NA , ND 

2 ,TRR,TRA,TAA,RTRR,RTRA,  C, CPA, CD 

X ,ALPD(2) ,RDU,RDV 

REALMS  U(MM,5) ,BA(MM,2) ,RR(1 ) 

EMA  U,BA,RR 

COMMON 

1 /COEFF/  C(MM,2,2) ,CPA(MM,2) ,CD<MM,2) 

3 /DATA/  U<MM,2) ,V(MM,2) ,ALP<MM,2) ,P(MM) ,R(MM) 

4 /DRAGl/  RH0<4) ,MUEF<2) , V18(2) ,NA,ND 
6 /FANDS/  F(MM,5> ,S(MM,5) 

8 /TAU/  TRR<MM,2) ,TRA(MM,2) ,TAA(MM,2) ,RTRR (MM,2) ,RTRA(MM,2) 

»V  CALL  DWMPY<W, 1 ,U, 1 ,F, 1 ,2*MM) 

*V  CALL  DWMPY(U(1 ,3) , 1 ,U, 1 ,F( 1 ,3) , 1 ,2«MM) 

*V  CALL  DWMOV(U, 1 ,F< 1 ,5) , 1 ,MM) 

DO  20  J=1,N2 
F< J,5)=W< J, 1 ) 

S( J,5)=0 . 

RDU=RR ( J)*<U< J, 1 )-U< J,2) ) 

RDV=RR( J)»(V( J, 1 )-V< J,2) ) 

DO  20  K=1 ,2 
KP2=K+2 

F(J,K)=U(J,K)»U( J,K) 

F( J,KP2)=W( J,KP2)*U( J,K ) 

S( J,K )=ALP< J,K )»< V( J,K )**2+CD( J,K )«RDU-C( J ,K , 1 )»TAA( J, 1 ) 

1 -C( J,K ,2)»TAA( J,2) ) 

20  S< J,KP2)=ALP ( J,K )»<-U< J,K )»V<  J,K )+CD( J ,K )»RDU+C( J ,K , 1 )*<  BA( J, 1 ) 
1 +TRA( J, 1 ) )+C( J,K ,2)»(BA( J,2)+TRA( J,2) ) ) 


RETURN 

END 


SUBROUTINE  JET < BA , RHO , BR , V , UJ , NG ) 
C 

C INJECTION  MOMENTUM  SOURCE 

C 


28 


0779 

0780 

0781 

0782 

0783 

0784 

0785 

0786 

0787 

0788 

0789 

0790 

0791 

0792 

0793 

0794 

0795 

0796 

0797 

0798 

0799 

0800 

0801 

0802 

0803 

0804 

0805 

0806 

0807 

0808 

0809 

081  0 

0811 

0812 

0813 

0814 

0815 

0816 

0817 

0818 

0819 

0820 

0821 

0822 

0823 

0824 

0825 

0826 

0827 

0828 

0829 

0830 

0831 

0832 

0833 

0834 

0835 

0836 

0837 

0838 


PARAMETER  (MM=101) 

REAL»8  BA (MM, 2) ,RHO< 1 ) , BR ( 1 ) , V < MM , 2 > , VJ 
EMA  BA,BR,V 

DO  10  J=1,NG 
IF<  BR(J)  .GT.  0.)  THEN 

Q=BR(J)»yj  ! Flow  rare/radian 

BA( J, 1 )=Q»(VJ-V( J, 1 ) )»RHO< 1 > ! NET  MOMENTUM  GAIN 

BA( J,2)=Q»<VJ-V( J,2) )*RH0<2) 

ELSE 

BA( J,1 )=0 .DO 
BA< J,2)=0 . DO 
ENDIF 

10  CONTINUE 
RETURN 
END 


*»**»«**»»»»*»»»*»*»***«*»***»»**»»»*»»**«»»»*»»*»**«»*****»**»■»**»**» 
SUBROUTINE  SIZES(D1 ,D2,ALP1 >,<860425. 1537) 

C TO  DETERMINE  THE  PARTICLE  DIAMETERS 

C 

REAL»8  D1 ,D2,ALP1 , DO, GAMMA 
COMMON  /DSIZE/  D0(2) ,GAMMA(2) 

Dl=DO(l >»ALP1**GAMMA(1 ) 

D2=D0(2)*(1 . ODO-ALPl )»*GAMMA(2) 

RETURN 

END 


«*«*««»««*««««««»««»««« »«««««»»»»*««»»*« »»««*«« »*»««««»*«* 
SEMA  /DATA/,/TAU/ 

SUBROUTINE  TAUOFU < MU , DR , RR , N2 ) , <360425. 1537) 

C STRESSES  AND  THEIR  DERIVATIVES 

C 

REAL*8  MU<2) ,DR,RR(1) 

EMA  RR 


PARAMETER  (MM=101) 

REAL»8  U,V,ALP,P,R 
2 ,TRR,TRA,TAA,RTRR,RTRA 

X ,TAMU,TMU 


COMMON 

3 /DATA/  U(MM,2) ,V(MM,2) ,ALP(MM,2) ,P(MM) ,R(MM) 

8 /TAU/  TRR (MM,2) , TRA(MM,2) ,TAA<MM,2) ,RTRR (MM,2) ,RTRA(MM, 2) 


DO  50  K=1 ,2 

CALL  DERI VI <U(1 ,K ) ,TRR ( 1 ,K ) ,DR ,N2) 

»V  CALL  DWMPY(ALP( 1 ,K) , 1 ,TRR ( 1 ,K ) , 1 ,TRR( 1 ,K) , 1 ,N2) 

*V  CALL  DWSMY(2 . »MU(K ) ,TRR ( 1 ,K ) , 1 ,TRR ( 1 ,K ) , 1 ,N2) 

»V  CALL  DWDIV(V<1 ,K) ,1 ,RR,1 ,TAA(1 ,K) ,1 ,N2) ! V/R 

DO  1 0 J=1 ,N2 

10  TAA( J,K )=V< J,K )/RR < J) 

CALL  DERI VI (TAA< 1 , K > , TR A ( 1 , K ) , DR , N2 ) 


I dU/dR 
! »ALP 
! »2*MU=TRR 


! dO/dR 


TMU=2 . *MU<K ) 

DO  20  J=1 ,N2 
TAMU=TMU*ALP( J,K) 

TRR ( J,K )=TAMU*TRR( J,K) 
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0839 

0840 

0841 

0842 

0843 

0844 

0845 

0846 

0847 

0848 

0849 

0850 

0851 

0 852 

0853 

0854 

0855 

0856 

0857 

0858 

0859 

0860 

0861 

0862 

0863 

0864 

0365 

0866 

0867 

0868 

0869 

0870 

0871 

0872 

0873 

0874 

0875 

0876 

0877 

0878 

0879 

0880 

0881 

0882 

0883 

0884 

0885 

0886 

0887 

0888 

0889 

0890 

0891 

0892 

0893 

0894 

0895 

0896 

0897 

0898 


TR A ( J , K ) =MU ( K ) »ALP ( J , K ) *RR ( J ) *TR A ( J , K ) 
TAA ( J , K ) =TAMU»U ( J , K ) /R  R ( J ) 

RTRR<J,K )=RR< J)*TRR< J,K ) 

20  RTRA<J,K)=RR(J)»TRA(J,K) 


*v 

»v 

*v 

CALL 

CALL 

CALL 

DWMPY(RR, 1 , TRAC  1 , K ) , 1 , TR A ( 1 , K ) , 1 , N2 ) 

DUMPY (ALP (1 ,K) ,1 ,TRA(1,K) ,1 , TRACI ,K),1 ,N2) 
DWSMYCMUCK) , TRACI , K ) , 1 , TR A ( 1 , K ) , 1 , N2 ) 

! «R 
! *ALP. 

' «MU=TRA 

»v 

*v 

»v 

CALL 

CALL 

CALL 

DUDIVCUC 1 ,K ) , 1 ,RR , 1 ,TAA( 1 ,K ) , 1 ,N2) ! U/R 
DWMPYCALPC 1 ,K) , 1 ,TAA< 1 ,K ) , 1 ,TAA( 1 ,K ) , 1 ,N2) 
DUSMYC2 . »MU(K ) ,TAA( 1 ,K ) , 1 ,TAA( 1 ,K ) , 1 ,N2) 

I «ALP 
! »2*MU=TAA 

»v 

»v 

50 

CALL  DWMPYCRR , 1 ,TRR(1 ,K ) , 1 ,RTRR( 1 ,K ) , 1 ,N2) 
CALL  DWMPYCRR , 1 , TR A < 1 , K ) , 1 , RTR A ( 1 , K ) , 1 , N2 ) 
CONTINUE 
RETURN 
END 

! R»TRR 
' R»TRA 

$EMA  /DATA/ 

SUBROUTINE  UOFU < W , RR , N2 ) ,<860425.1537) 

C CONVERTS  U TO  THE  INDEPENDENT  VAR  TABLES ( U , V , ALP ) 

PARAMETER  (MM=101) 

REAL»8  U(MM,5)*RR(MM) 

EMA  W ,RR . 

REAL*8  ALMT,  U,V,ALP,P,R 
X ,UJ6 

COMMON 

X /ALPLMT/  ALMT(2) 

3 /DATA/  U(MM,2) , V(MM,2) , ALP (MM,2> ,P <MM)  ,R  (MM) 

»V  CALL  DUDIV(U( 1 ,5) , 1 ,RR , 1 , ALP , 1 ,N2) 

C CHECK  VOLUME  FRACTION  6 FLOW  DIRECTIONS 

DO  50  J=1 ,N2 
ALP( J, 1 )=W( J,5)/RR  < J) 

IF(ALP<J,1)  .LE.  ALMT(l)  .OR.  ALP(J,1)  . GE . ALMT(2>)  THEN 
IF(ALP(J,1)  .LE.  ALMT(l))  THEN 
ALP( J, 1 )=ALMT< 1 ) 

W< J,3)=U( J,4)*ALMT(1 )/( 1 . D0-ALMT< 1 ) ) 

ELSE 

ALP( J, 1 )=ALMT(2) 

W<J,4)=U(J,3)»(1 .D0-ALMT<2) ) /ALMT (2) 

ENDIF 

W( J, 1 )»0 .DO 
W<J,2)»0 .DO 
W< J,5)=ALP( J, 1 )*RR< J) 

ENDIF 

ALP( J,2)  = l .D0-ALP( J,  1 ) 

IF(  W(J,2)  .LT.  0.)  THEN  ' PHASE-2  DOES  NOT  MOVE  IN 

W< J,2)=0 . 

W( J, 1 ) = 0 . 

ENDIF 
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0899 

0900 

0901 

0902 

0903 

0904 

0905 

0906 

0907 

0908 

0909 

0910 

0911 

0912 

0913 

0914 


U< J,1 )=W<J,1 )/U<J,5) 
V( J, 1 )=W< J,3)/U< J,5) 
WJ6=RR( J)»ALP<J,2) 

U( J,2)=W( J,2)/WJ6 
J,2)=W< J,4)/UJ6 
CONTINUE 


CALL  DUDIU(W,1 ,U( 1 ,5) , 1 ,U, 1 ,N2) 

CALL  DUDIV<W(1 , 3 ) , 1 , W (1 , 5 ) , 1 , V , 1 , N2 ) 
CALL  DWHPY(ALP( 1 ,2),1,RR,1,U<1,6),1, 
CALL  DWDIU(U(1 ,2) ,1 ,U<1,6),1,U(1,2), 
CALL  DUDIV(W<1 ,4),1,U(1,6),1,V(1,2), 

RETURN 

END 


N2) 

1 ,N2) 
1 ,N2) 
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Exhibit  A 


:GLVM 

Enter 

90 

A Sample  Input 

lu.  for  saving  data.  D . F . = b 

Enter  hILE  NAME  tor  saving  data. 
TS153r :LB 

Enter  N0TES<<73  CHAR.)  for  the  job 
SAMPLE  RUN  OF  TEST  153 


enter 

Reynolds  no., RE.  D.F.=  1.0000E+05 

« 

Enter 

D,F.= 

DENSITY  and  VISCOSITY  ratios.  ' . 

1 .2930E-03  1 . 1295E-02 

Enter 

D.F.= 

BASE  DIAMETERS:  D01 ,D02 
1 .OOOOE-02  1 .OOOOE-02 

f * 

Enter 
D.  F. 

SIZE  EXPONENT:  GAMMA  1 ,GAMMA2 
= 2.0000E-01  2.0000E-01 

* t 

Enter 

ueighting  exponent:  NA,ND.  D.F.=  *4  . OOOOE  t-OO  ^.'JOOOE  + OU 

« « 

Enter 

D.F.= 

GAS  VOLUME  FRACTION  1 imi ts : ALMT1  . ALMT2  . 
1 .OOOOE-04  9.9999E-01 

Enter  eddy  visicosity  factor.  D.F.=  1 . 1 4'^8E■'-‘J2  I.OuQOEt-03 

1000,1000 


Enter 

wall  condi  t ion  , i.  = nonsl  ip  . 0 = si  IP  . D.t.  ’ 

Enter 

nujTiericai  damping  factor.  D . F . =■•  0.0000E-*-00 

Enter 

JET  SIZE  defined  by  RJ1.RJ2.  D.F.=  3.5000E-01  9.5000E-0 

« * 

Enter 
D.F,= 
1 ,0,1 
Enter 
D.F.= 

INJECTION  SPEED  AND  TIME  RANGE.  VJ.T1,T2 
1.0000E+01  O.OOOOE+00  1.0000E+01 

data  FIl.E  iNAME  for  initial  cond  . i^  any: 
Simple  vortex 

Enter 

initial  value  of  alpl.  D.F.=  2.5000E-t.M 

Enter  type  of  vortex:  0=At  rest.1=pure  rotation 
2=H.O. ,3=GIT.  D.F.=  0 


Enter 

0.^01 

Enter 

.01 

1 NP  = 

1 NP  = 

INITIAL  and  FINAL  TIMEs.  D.F.=  O.UOOOE^OO  5.00UOE+00 
TIME  STEP  for  output.  D.F.=  2.0000E-01 

1 T=0.00GE+00  NT-  2 DT-i.OOOE-06 

2 T = l.200E-02  NT=  b DT-3.000F.-n3 
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Exhibit  B 


A Sample  Output 


TS153  T--00004  13  ON  CR  LB  USING  00024  PLKS  R = 0000 

0001  1 SAMPLE  RUN  OF  TEST  153 

0002  INITIAL  CONDITION  FILE:  . 


0003 

^♦DIMENSION  UNITS  ARE  IN 

MKS** 

0004 

DENSITY  SCALE(kg/fn3) 

! . OOnOE  + 03 

0005 

LENGTH  SCALE=RTANK,(m) 

1 .OOOOE+00 

0006 

VELOCITY  SCALE<m/s) 

1 .5I40E-01 

0007 

TIME  SACLE(s) 

6.6050E+00 

0008 

0009 

PRESSURE  SCALE(Pa) 

2.2922E>01 

0010 

Reynolds  number,  Re 

1 .0000E^05 

001  1 

Jet  Size,  RJ1 , RJ2 

8.5000E-01 

9.5000E-01 

0012 

Tangential  jet,  Q J , V J 

1 .OOOOE-01 

1 .OOOOE+00 

0013 

0014 

Injection  time,TJ1,TJ2 

O.OOOOE+00 

1 .OOOOE+00 

0015 

PHASE- 1 

PHASE-2 

0016 

Dens  1 ty 

1 .2930E-03 

1 . OOOOE  + OO 

0017 

U iscos  i ty 

1 . 1295E-07 

1 .OOOOE-05 

0018 

Eddy  viscosity  factor 

1 .OOOOE+03 

1 .0000E  + 03 

0019 

Base  dia. 

1 .OOOOE-02 

1 .OOOOE-02 

0020 

S 1 ze  exp . 

2.0000E-01 

2.0000E-01 

0021 

Phase  limits 

1 .OOOOE-04 

T .0014E-05 

0022 

0023 

0024  OTHER  CONSTANS:  IN , IVTX . NA , ND . DAMP . VPEAK , RPEAK , OMEGA , D 1 /D2 . MU  I /Mi 


0025 

1 0 

4. 

4.  0.00 

0026 

0. OOOOE+OO 

1 .OOOOE-01 

0. OOOOE+OO 

1 .2930E-03 

1 . 1295E-02 

0027 

0028 

101 

. 10 

0029 

1 NP  = 

1 

T=0.000E+00  NT=  2 

DT=1  .OOOE- 

06 

0030 

J 

ALPi 

U1 

U2 

VI 

V2 

0031 

1 

.2500 

O.OOOE+00 

O.OOOE^OO 

O.OOOE+00 

O.OOOE+00 

0 .OOOE  + 0 

0032 

2 

.2500 

O.OOOE+00 

O.OOOE+00 

O.OOOE+00 

O.OOOE+00 

o.oooE+n 

QUTP 

UT  IN  THE  BETWEEN  OMITTED  

0230 

97 

.2500 

-8.685E- 10 

2.895E-10 

3.513E-04 

3.474E-04 

3.969E-0 

0231 

98 

.2500 

-1 .775E-10 

5.918E- 1 I 

1 .268E-04 

1 .250E-04 

3.972E-0 

0232 

99 

.2500 

-2.785E-1 1 

9.284E-12 

4.319E-05 

4.256E-05 

3.973E-0 

0233 

100 

.2500 

-5.225E-12 

1 .742E-12 

1 .553E-05 

1 .531E-05 

3.373E-0 

0234 

101 

.2500 

O.OOOE+00 

O.OOOE+00 

O.OOOE+00 

O.OOOE+00 

3.973E-n 

THE  REST  OF  THE  OUTPUT  IS  OMITTED 
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Complete  etep 
Correction  step 


Prediction  step 


Initial  value 


Figure  1.  Two  Step  Difference  Scheme  (Backward  Predictor  - 
Forward  Corrector  Version) 
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Figure  2.  Velocity  Vector  Distributions 

The  annular  region  between  two  dashed  circles  Is 
the  region  of  Injection. 
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